      subroutine vinit
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use boundary
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M, ONLY:
      use cophys_com_M
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx
      use blcom_com_M
      use objects_com_M, ONLY: nrg_obj, chi, chic
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer :: n, ijkb, ijkt, ijkl, ijkr, is
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double) :: phibc, qtotl, qerror, rvol, wk, dummy, dto2, rvvol, &
         divtotl, qvtot, wsp, rcntr, turnon, jtmp
      real(double) :: tiny, error, relax
      logical :: expnd, nomore
      character :: name*80
!-----------------------------------------------
!
      phibc = 0.
!
!     a routine to prepare varibles for the next computationcycle
!
!     normalize the pressure
!
      qtotl = 0.0
      qerror = 0.0
      do n = 1, ncells
         ijk = ijkcell(n)
         rvol = 1./vol(ijk)
         qdnc(ijk) = qdnc(ijk)*rvol
!
         pixx(ijk) = pixx(ijk)*rvol
         pixy(ijk) = pixy(ijk)*rvol
         pixz(ijk) = pixz(ijk)*rvol
         piyy(ijk) = piyy(ijk)*rvol
         piyz(ijk) = piyz(ijk)*rvol
         pizz(ijk) = pizz(ijk)*rvol
!
         qtotl = qtotl + qdnc(ijk)*vol(ijk)
!
      end do
      write (*, *) 'vinit: qtotl from parcel', qtotl

      if (nrg_obj > 0) then

         do n = 1, ncells
            ijk = ijkcell(n)
            rvol = 1./vol(ijk)
            chic(ijk) = chic(ijk)*rvol
         end do

      endif
!
!      call bccz(ibar+2,jbar+2,kbar+2,1.,pixx)
!      call bccz(ibar+2,jbar+2,kbar+2,1.,piyy)
!      call bccz(ibar+2,jbar+2,kbar+2,1.,pizz)
!      call bccz(ibar+2,jbar+2,kbar+2,1.,pixy)
!     the following should have -1
!      call bccz(ibar+2,jbar+2,kbar+2,-1.,pixz)
!      call bccz(ibar+2,jbar+2,kbar+2,-1.,piyz)
!
      wk = 1.
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, qdnc, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, pixx, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, pixy, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, pixz, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, piyy, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, piyz, phibc, wk, a11, periodicx)
      call ctovtoc (nvtx, ijkvtx, ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1&
         , kbp1, nsmooth, pizz, phibc, wk, a11, periodicx)
!
!      call torusbc(ibar+2,jbar+2,kbar+2,cdlt,sdlt,
!     &     dummy,dz,pixx,piyy,pizz,periodicx)
!      call torusbc(ibar+2,jbar+2,kbar+2,cdlt,sdlt,
!     &     dummy,dz,pixy,pixz,piyz,periodicx)
!
!      call bccz(ibar+2,jbar+2,kbar+2,1.,pixx)
!      call bccz(ibar+2,jbar+2,kbar+2,1.,piyy)
!      call bccz(ibar+2,jbar+2,kbar+2,1.,pizz)
!      call bccz(ibar+2,jbar+2,kbar+2,1.,pixy)
!     the following should have -1
!      call bccz(ibar+2,jbar+2,kbar+2,-1.,pixz)
!      call bccz(ibar+2,jbar+2,kbar+2,-1.,piyz)
!
      call divpi (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, pixx, pixy, pixz, piyy, piyz, pizz, divpix, &
         divpiy, divpiz)
!
!     impose toroidal boundary conditions on vertex variables
!
      dummy = 0.0
!
      call torusbcv (ibp1 + 1, jbp1 + 1, kbp1 + 1, cdlt, sdlt, dummy, dz, &
         jxtilde, jytilde, jztilde, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, jxtilde, jytilde, jztilde)
!
      call torusbcv (ibp1 + 1, jbp1 + 1, kbp1 + 1, cdlt, sdlt, dummy, dz, &
         j12x, j12y, j12z, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, j12x, j12y, j12z)
!
      call torusbcv (ibp1 + 1, jbp1 + 1, kbp1 + 1, cdlt, sdlt, dummy, dz, &
         divpix, divpiy, divpiz, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, divpix, divpiy, divpiz)
!
!     impose reflection bc in z
!
!      if(cartesian) then
!
      dummy = 2.
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, jxtilde)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, jytilde)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, jztilde)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, j12x)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, j12y)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, j12z)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, divpix)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, divpiy)
!
      call bcvz (ibar + 2, jbar + 2, kbar + 2, dummy, divpiz)
!
!      end if
!
      do k = 2, kbp2
!
         i = 2
         if (ibp2 - 1 > 0) then
            do i = 1, ibp2 - 1
!
               jmag(iwid*i+1+jwid+(k-1)*kwid) = jmag(iwid*i+1+jwid+(k-1)*kwid)&
                   + jmag(iwid*i+1+jbp1*jwid+(k-1)*kwid)
               jmag(iwid*i+1+jbp1*jwid+(k-1)*kwid) = jmag(iwid*i+1+jwid+(k-1)*&
                  kwid)
!
            end do
            i = ibp2 + 1
         endif
!
         if (periodicx) then
!
            j = 2
            if (jbp2 - 1 > 0) then
               do j = 1, jbp2 - 1
!
                  jmag(jwid*j+1+iwid+(k-1)*kwid) = jmag(jwid*j+1+iwid+(k-1)*&
                     kwid) + jmag(jwid*j+1+ibp1*iwid+(k-1)*kwid)
                  jmag(jwid*j+1+ibp1*iwid+(k-1)*kwid) = jmag(jwid*j+1+iwid+(k-1&
                     )*kwid)
!
               end do
               j = jbp2 + 1
            endif
!
         else
!
            j = 2
            if (jbp2 - 1 > 0) then
               j = jbp2 + 1
            endif
!
         endif
!
      end do
!
!
!     include pressure contribution to the current
!     and magnetization current
!
      dto2 = 0.5*dt
      do n = 1, nvtx
!c
         ijk = ijkvtx(n)
!c
         rvvol = 1./vvol(ijk)
         jmag(ijk) = jmag(ijk)*rvvol
!c
         jxtilde(ijk) = jxtilde(ijk)*rvvol
         jytilde(ijk) = jytilde(ijk)*rvvol
         jztilde(ijk) = jztilde(ijk)*rvvol
!c
         j12x(ijk) = j12x(ijk)*rvvol
         j12y(ijk) = j12y(ijk)*rvvol
         j12z(ijk) = j12z(ijk)*rvvol
!c
!         pixx(ijk) = jxtilde(ijk)
!         pixy(ijk) = jytilde(ijk)
!         pixz(ijk) = jztilde(ijk)
!c
!      jxtilde(ijk)=jxtilde(ijk)-(divpix(ijk)*dto2
!     &     +jmag(ijk)*gradxb(ijk)*dto2)*RVVOL
!      jytilde(ijk)=jytilde(ijk)-(divpiy(ijk)*dto2
!     &     +jmag(ijk)*gradyb(ijk)*dto2)*RVVOL
!      jztilde(ijk)=jztilde(ijk)-(divpiz(ijk)*dto2
!     &     +jmag(ijk)*gradzb(ijk)*dto2)*RVVOL
!c
         jxtilde(ijk) = jxtilde(ijk) - divpix(ijk)*dto2*rvvol
         jytilde(ijk) = jytilde(ijk) - divpiy(ijk)*dto2*rvvol
         jztilde(ijk) = jztilde(ijk) - divpiz(ijk)*dto2*rvvol
!c
         jtmp=j12x(ijk)
         j12x(ijk)=.5d0*(j12x(ijk)+joldx(ijk))
         joldx(ijk)=jtmp
         jtmp=j12y(ijk)
         j12y(ijk)=.5d0*(j12y(ijk)+joldy(ijk))
         joldy(ijk)=jtmp
         jtmp=j12z(ijk)
         j12z(ijk)=.5d0*(j12z(ijk)+joldz(ijk))
         joldx(ijk)=jtmp
!c
      end do
!
!
!     Boundary Cond. in z,x
!
!      call bcvz(ibar+2,jbar+2,kbar+2,0.,jztilde)
!
!      call torusbc_scal(2,ibp2,2,jbp2,2,kbp2,iwid,jwid,kwid,
!     &     0.,jxtilde,periodicx)
!      call torusbc_scal(2,ibp2,2,jbp2,2,kbp2,iwid,jwid,kwid,
!     &     0.,jztilde,periodicx)
!
      wk = 1.
      call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, nsmooth, &
         sdlt, cdlt, wkl, wkr, wkf, wke, phibc, jxtilde, jytilde, jztilde, a11&
         , a12, a13, periodicx)
      wk = 1.
      call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, nsmooth, &
         sdlt, cdlt, wkl, wkr, wkf, wke, phibc, j12x, j12y, j12z, a11&
         , a12, a13, periodicx)
!
!     Boundary Cond. in z,x
!
!      call bcvz(ibar+2,jbar+2,kbar+2,0.,jztilde)
!
!      call torusbc_scal(2,ibp2,2,jbp2,2,kbp2,iwid,jwid,kwid,
!     &     0.,jxtilde,periodicx)
!      call torusbc_scal(2,ibp2,2,jbp2,2,kbp2,iwid,jwid,kwid,
!     &     0.,jztilde,periodicx)
!

!
!     if implicit field solve, calculate estimate of charge at dto2
!
      if (.not.explicit) then
!
         call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
            , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
            c3z, c4z, c5z, c6z, c7z, c8z, vol, jxtilde, jytilde, jztilde, piyy)
!
         divtotl = 0.0
         do n = 1, ncells
!
            ijk = ijkcell(n)
!
            qdnc0(ijk) = qdnc(ijk)
            qdnc(ijk) = qdnc(ijk) - piyy(ijk)*cntr*dt
!
            divtotl = divtotl + piyy(ijk)*vol(ijk)*cntr*dt
!      divtotl=divtotl+qdnc(ijk)*vol(ijk)
!
         end do
!

         write (*, *) 'vinit: divtotal ', divtotl
!
      endif
!
!     Set bc on qdnv in z
!
      qvtot = 0.
      do is = 1, nsp
!
         wsp = 1.
         call torusbc_scal (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, wsp, &
            qdnv(1,is), periodicx)
!
         if (cartesian) call bcvz (ibar + 2, jbar + 2, kbar + 2, 2., qdnv(1,is)&
            )
!
!
         do n = 1, nvtx
            ijk = ijkvtx(n)
            rvvol = 1./vvol(ijk)
            qvtot = qvtot + qdnv(ijk,is)
            qdnv(ijk,is) = qdnv(ijk,is)*rvvol
         end do
      end do

      write (*, *) 'qvtot:', qvtot
!
      a21(:ibp2*jbp2*kbp2) = 0.
      a22(:ibp2*jbp2*kbp2) = 0.
      do is = 1, nsp
         wk = 1.
         call vtoctov (nvtx, ijkvtx, iwid, jwid, kwid, ibp2, jbp2, kbp2, &
            nsmooth, sdlt, cdlt, wkl, wkr, wkf, wke, phibc, qdnv(1,is), a21, &
            a22, a11, a12, a13, periodicx)
      end do

!
!      write (*, *) 'vinit:renormalize udrift'
!      do n = 1, nvtx
!         ijk = ijkvtx(n)
!
!         rvvol = 1./vvol(ijk)
!
!         pixx(ijk) = pixx(ijk)*rvvol
!         pixy(ijk) = pixy(ijk)*rvvol
!         pixz(ijk) = pixz(ijk)*rvvol
!     treat susceptibility similarly to vertex charge density
!
!      end do
!
      if (nrg_obj > 0) then

         do n = 1, nvtx
            ijk = ijkvtx(n)
            rvvol = 1./vvol(ijk)
            chi(ijk) = 1. + chi(ijk)*rvvol
         end do

      endif

!
!     avergaing for graphics
!
      if(ncyc.le.2) then

      uxsp1=jxs(:,1)*(qdnv(:,1)+1.E-10)
      uysp1=jys(:,1)*(qdnv(:,1)+1.E-10)
      uzsp1=jzs(:,1)*(qdnv(:,1)+1.E-10)
      densp1(:)=qdnv(:,1)
      uxsp2=jxs(:,2)*(qdnv(:,2)+1.E-10)
      uysp2=jys(:,2)*(qdnv(:,2)+1.E-10)
      uzsp2=jzs(:,2)*(qdnv(:,2)+1.E-10)
      densp2(:)=qdnv(:,2)


      else

      uxsp1=uxsp1*avg+(1d0-avg)*jxs(:,1)*(qdnv(:,1)+1.E-10)
      uysp1=uysp1*avg+(1d0-avg)*jys(:,1)*(qdnv(:,1)+1.E-10)
      uzsp1=uzsp1*avg+(1d0-avg)*jzs(:,1)*(qdnv(:,1)+1.E-10)
      densp1=densp1*avg+(1d0-avg)*qdnv(:,1)
      uxsp2=uxsp2*avg+(1d0-avg)*jxs(:,2)*(qdnv(:,2)+1.E-10)
      uysp2=uysp2*avg+(1d0-avg)*jys(:,2)*(qdnv(:,2)+1.E-10)
      uzsp2=uzsp2*avg+(1d0-avg)*jzs(:,2)*(qdnv(:,2)+1.E-10)
      densp2=densp2*avg+(1d0-avg)*qdnv(:,2)


      exavp(:ibp2*jbp2*kbp2) = exavp(:ibp2*jbp2*kbp2)&
           *avg+ex(:ibp2*jbp2*kbp2)*(1d0-avg)
      eyavp(:ibp2*jbp2*kbp2) = eyavp(:ibp2*jbp2*kbp2)&
           *avg+ey(:ibp2*jbp2*kbp2)*(1d0-avg)
      ezavp(:ibp2*jbp2*kbp2) = ezavp(:ibp2*jbp2*kbp2)&
           *avg+ez(:ibp2*jbp2*kbp2)*(1d0-avg)

      end if
!
      return
      end subroutine vinit
